plus sn10x

load dependancies

Baf53cre ENS neurons, SI data
Nb infection 5d, PBS(control) x4 INF(inflammation) x4

loading 8k cells for each,
demo called 20,985 cells
plus called 21,294 cells

cleaning-up and re-clustering and annotation
ref-mapping to NatNeur2021

GEX

mainly follow https://satijalab.org/seurat/articles/pbmc3k_tutorial.html

load obj

GEX.seur <- readRDS("./GEX.seur.preAnno0717.rds")
GEX.seur

## An object of class Seurat 
## 23258 features across 9423 samples within 1 assay 
## Active assay: RNA (23258 features, 1500 variable features)
##  3 dimensional reductions calculated: pca, tsne, umap

color.FB <- ggsci::pal_d3("category20c")(20)[c(2,7,12,17,
                                               1,6,11,16
                                               )]

filtering

# Visualize QC metrics as a violin plot
VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.rb"), ncol = 4, pt.size = 0.01, group.by = "FB.info", cols = color.FB)

VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.rb"), ncol = 4, pt.size = 0, group.by = "FB.info", cols = color.FB)

plota <- FeatureScatter(GEX.seur, feature1 = "nCount_RNA", feature2 = "percent.mt", group.by = "FB.info", cols = color.FB) 
plotb <- FeatureScatter(GEX.seur, feature1 = "nCount_RNA", feature2 = "nFeature_RNA", group.by = "FB.info", cols = color.FB) 
plotc <- FeatureScatter(GEX.seur, feature1 = "nCount_RNA", feature2 = "percent.rb", group.by = "FB.info", cols = color.FB)
plota + plotb + plotc

subset

head(GEX.seur@meta.data)

##                       orig.ident nCount_RNA nFeature_RNA percent.mt percent.rb
## AAACCCACAAGACGAC-1 Baf53Nb_Ileum       3257         1801 0.36843721  0.3991403
## AAACCCAGTACAGGTG-1 Baf53Nb_Ileum       5446         2500 0.09181050  0.2937936
## AAACCCAGTGGGCTCT-1 Baf53Nb_Ileum       1511          997 0.66181337  0.4632694
## AAACCCAGTTTGTTCT-1 Baf53Nb_Ileum       2855         1577 0.98073555  0.3152364
## AAACCCATCCTAGCCT-1 Baf53Nb_Ileum       2433         1451 0.08220304  0.3699137
## AAACCCATCGAAACAA-1 Baf53Nb_Ileum       3129         1656 0.12783637  0.4474273
##                        FB.info      S.Score     G2M.Score Phase RNA_snn_res.1.5
## AAACCCACAAGACGAC-1 Baf53Nb.IF4  0.011590405 -0.0004169865     S               3
## AAACCCAGTACAGGTG-1 Baf53Nb.CL3 -0.038183546  0.0022744719   G2M               0
## AAACCCAGTGGGCTCT-1 Baf53Nb.CL4 -0.024203070  0.0012414826   G2M               7
## AAACCCAGTTTGTTCT-1 Baf53Nb.IF1 -0.013980476  0.0039329410   G2M               3
## AAACCCATCCTAGCCT-1 Baf53Nb.IF2 -0.028925620 -0.0132582758    G1               0
## AAACCCATCGAAACAA-1 Baf53Nb.CL3 -0.008077289 -0.0028336129    G1              13
##                    seurat_clusters sort_clusters preAnno1 preAnno2        cnt
## AAACCCACAAGACGAC-1               3             3   PEMN.5     PEMN Baf53Nb.IF
## AAACCCAGTACAGGTG-1               0             0   PEMN.2     PEMN Baf53Nb.CL
## AAACCCAGTGGGCTCT-1               7             7    PSN.1      PSN Baf53Nb.CL
## AAACCCAGTTTGTTCT-1               3             3   PEMN.5     PEMN Baf53Nb.IF
## AAACCCATCCTAGCCT-1               0             0   PEMN.2     PEMN Baf53Nb.IF
## AAACCCATCGAAACAA-1              13            13    PIN.2      PIN Baf53Nb.CL
##                    pANN_0.25_0.005_471 DoubletFinder0.05 pANN_0.25_0.005_942
## AAACCCACAAGACGAC-1          0.07936508           Singlet          0.07936508
## AAACCCAGTACAGGTG-1          0.58730159           Doublet          0.58730159
## AAACCCAGTGGGCTCT-1          0.00000000           Singlet          0.00000000
## AAACCCAGTTTGTTCT-1          0.03174603           Singlet          0.03174603
## AAACCCATCCTAGCCT-1          0.01587302           Singlet          0.01587302
## AAACCCATCGAAACAA-1          0.00000000           Singlet          0.00000000
##                    DoubletFinder0.1
## AAACCCACAAGACGAC-1          Singlet
## AAACCCAGTACAGGTG-1          Doublet
## AAACCCAGTGGGCTCT-1          Singlet
## AAACCCAGTTTGTTCT-1          Singlet
## AAACCCATCCTAGCCT-1          Singlet
## AAACCCATCGAAACAA-1          Singlet

table(GEX.seur@meta.data[,c("preAnno1","DoubletFinder0.05")])

##         DoubletFinder0.05
## preAnno1 Doublet Singlet
##   PEMN.1       0     799
##   PEMN.2      22    1076
##   PEMN.3      22     697
##   PEMN.4      27     458
##   PEMN.5       5     739
##   PEMN.6      16     287
##   PIMN.1       4     923
##   PIMN.2       6     567
##   PIMN.3       6     561
##   PIMN.4       8     223
##   PIN.1       10     182
##   PIN.2       25     243
##   PIN.3        7     100
##   PSN.1       36     457
##   PSN.2       33     335
##   PSN.3       20     204
##   PSN.4       38     409
##   PSN.5        7      40
##   PSVN.1       1      39
##   PSVN.2      11     146
##   PSVN.3      26     256
##   MIX.1      121      50
##   MIX.2       14      99
##   Glia         6      62

table(GEX.seur@meta.data[,c("preAnno1","DoubletFinder0.1")])

##         DoubletFinder0.1
## preAnno1 Doublet Singlet
##   PEMN.1       1     798
##   PEMN.2      61    1037
##   PEMN.3      56     663
##   PEMN.4      64     421
##   PEMN.5      13     731
##   PEMN.6      44     259
##   PIMN.1      14     913
##   PIMN.2      24     549
##   PIMN.3      27     540
##   PIMN.4      36     195
##   PIN.1       24     168
##   PIN.2       47     221
##   PIN.3       29      78
##   PSN.1       52     441
##   PSN.2       51     317
##   PSN.3       31     193
##   PSN.4       56     391
##   PSN.5       20      27
##   PSVN.1      16      24
##   PSVN.2      25     132
##   PSVN.3      38     244
##   MIX.1      138      33
##   MIX.2       42      71
##   Glia        33      35

GEX.seur <- subset(GEX.seur, subset = DoubletFinder0.1=="Singlet" & preAnno2 %in% c("PEMN","PIMN","PIN","PSN","PSVN"))
GEX.seur

## An object of class Seurat 
## 23258 features across 8342 samples within 1 assay 
## Active assay: RNA (23258 features, 1500 variable features)
##  3 dimensional reductions calculated: pca, tsne, umap

DimPlot(GEX.seur, reduction = "umap", group.by = "preAnno1", label = T, label.size = 2.5) + 
  DimPlot(GEX.seur, reduction = "umap", group.by = "preAnno2", label = T, repel = F)

FeaturePlot(GEX.seur, reduction = "umap", features = c("nFeature_RNA","nCount_RNA","percent.mt","percent.rb"))

# Visualize QC metrics as a violin plot
VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.rb"), ncol = 4, pt.size = 0.01, group.by = "FB.info", cols = color.FB)

VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.rb"), ncol = 4, pt.size = 0, group.by = "FB.info", cols = color.FB)

plota <- FeatureScatter(GEX.seur, feature1 = "nCount_RNA", feature2 = "percent.mt", group.by = "FB.info", cols = color.FB) 
plotb <- FeatureScatter(GEX.seur, feature1 = "nCount_RNA", feature2 = "nFeature_RNA", group.by = "FB.info", cols = color.FB) 
plotc <- FeatureScatter(GEX.seur, feature1 = "nCount_RNA", feature2 = "percent.rb", group.by = "FB.info", cols = color.FB)
plota + plotb + plotc

par(mar=c(6,4,2,3))
sl_stat <- table(GEX.seur$FB.info)
barplot(sl_stat,ylim = c(0,2300),
        col = c("#FF6C91","lightgrey",color.FB),
        main = "Feature Barcode statistics",cex.names = 0.75, xaxt = "n")
axis(1,1:10*1.2-0.48 ,levels(GEX.seur$FB.info), las=3, cex.axis=0.85)
text(x=1:10*1.2-0.45,y=sl_stat+125,paste0(sl_stat,"\n",100*round(as.numeric(sl_stat/sum(sl_stat)),4),"%"),cex = 0.75)

re-clustering

Normalizing

#GEX.seur <- NormalizeData(GEX.seur, normalization.method = "LogNormalize", scale.factor = 10000)

GEX.seur <- FindVariableFeatures(GEX.seur, selection.method = "vst", nfeatures = 1500)

# Identify the 10 most highly variable genes
#top10 <- head(VariableFeatures(GEX.seur), 10)
top20 <- head(VariableFeatures(GEX.seur), 20)

# plot variable features with and without labels
plot1 <- VariableFeaturePlot(GEX.seur)
plot2 <- LabelPoints(plot = plot1, points = top20, repel = T)
plot1 + plot2

head(VariableFeatures(GEX.seur), 100)

##   [1] "Mgat4c"        "Gal"           "Zfp804a"       "Grm5"         
##   [5] "Cntn5"         "Nmu"           "Klhl1"         "Adgrg6"       
##   [9] "Gm32647"       "4930432L08Rik" "Col24a1"       "Cpne4"        
##  [13] "Kcnb2"         "Gm30613"       "Robo2"         "Ntng1"        
##  [17] "Cntnap2"       "Tmeff2"        "Brinp3"        "Gm29521"      
##  [21] "Nrxn3"         "Lingo2"        "Ano2"          "Gpr149"       
##  [25] "Dgkg"          "Cdh8"          "Cdh18"         "Zfp804b"      
##  [29] "Gm21847"       "Ebf1"          "Cdh9"          "Ntrk3"        
##  [33] "Gm49953"       "Gm15261"       "Dgki"          "Sgcz"         
##  [37] "Pdzrn4"        "Cmah"          "Nxph2"         "Unc5d"        
##  [41] "Astn2"         "Kcnip4"        "1700111E14Rik" "Sema5a"       
##  [45] "Pcdh11x"       "Schip1"        "Pcdh9"         "Prkg1"        
##  [49] "Itgb6"         "Vwc2l"         "Cdh6"          "Nrg1"         
##  [53] "Postn"         "Dcc"           "Calcb"         "4930509J09Rik"
##  [57] "Car10"         "Piezo2"        "Grin3a"        "Myl1"         
##  [61] "March1"        "Clstn2"        "Vip"           "Gm38405"      
##  [65] "4930422I22Rik" "Rasgef1b"      "Bmpr1b"        "Efr3a"        
##  [69] "Ccbe1"         "Rab27b"        "Plxna4"        "Fut9"         
##  [73] "Dpp10"         "Gpc5"          "Hbb-y"         "AC150683.1"   
##  [77] "Gm49938"       "Gm20754"       "Gm30094"       "Pbx3"         
##  [81] "Gm15680"       "Asic2"         "Egfem1"        "9530059O14Rik"
##  [85] "Pcdh10"        "Abi3bp"        "Cysltr2"       "Gm15584"      
##  [89] "Serpini1"      "Col25a1"       "Skap1"         "Nkain3"       
##  [93] "Robo1"         "Fgf14"         "Gm28494"       "Trhde"        
##  [97] "Tafa2"         "Spock3"        "Chrna9"        "Zfhx3"

GEX.seur <- ScaleData(GEX.seur, features = rownames(GEX.seur))

## Centering and scaling data matrix

PCA

# exclude MT genes  and more 

#DIG <- grep("^Tra|^Trb|^Trg|^Trd|^Tcr|^Igm|^Igh|^Igk|^Igl|Jchain|^Hsp|^Rps|^Rpl|Hbb-|Hba-|^Dnaj|^AY|^Gm|^Hist",rownames(GEX.seur),value = T)

DIG <- grep("^Tra|^Trb|^Trg|^Trd|^Tcr|^Igm|^Igh|^Igk|^Igl|Jchain|Mzb1|Vpreb|Lars2|Jun|Fos|^Hsp|^Rps|^Rpl|Hbb-|Hba-|^Dnaj|^AY|^Gm|^Hist|Rik$|^AC|-ps",
            rownames(GEX.seur),value = T)
CC_gene <- Hmisc::capitalize(tolower(as.vector(unlist(cc.genes.updated.2019))))
MT_gene <- grep("^mt-",rownames(GEX.seur),value = T)

GEX.seur <- RunPCA(GEX.seur, features = setdiff(VariableFeatures(object = GEX.seur),
                                                c(MT_gene,
                                                  DIG, 
                                                  CC_gene) ))

## PC_ 1 
## Positive:  Ntrk3, Tmeff2, Robo2, Cdh8, Ano2, Nrxn3, Cpne4, Clstn2, Myl1, Mgat4c 
##     Dgkg, Gpr149, Pdzrn4, Adgrg6, Plxna4, Grin3a, Pcdh10, Cdh6, Spock3, Itgb6 
##     Cysltr2, Cntn5, Zfp804a, Sgcz, Iqgap2, Cacna2d3, Cntnap2, Ccbe1, Astn2, Cux2 
## Negative:  Lrrtm4, Galntl6, Tox, Ndst4, Cacnb2, Grik1, Kcnd2, Lrp1b, Bnc2, Pcdh7 
##     Epha5, Synpo2, Kcnab1, Syt6, Tshz2, Zfp536, Plcxd3, Cdc14a, Man2a1, Chrna7 
##     Specc1, Lrrc4c, Nos1, Galnt17, Pitpnc1, Lrrc7, Tenm3, Fat3, Dach1, Brinp2 
## PC_ 2 
## Positive:  Rbfox1, Bnc2, Ptprt, Tafa1, Tshz2, Grik1, Gpc6, Frmd4b, Plcxd3, St6galnac3 
##     Tox, Caln1, Brinp2, Fbxw15, Dock2, Cdc14a, Tcf7l2, Fbxw24, Tmem132c, Pcdh7 
##     Pld5, Oprk1, Sdk2, Specc1, Tafa5, Adamtsl1, Slc35f4, Nfia, Unc5c, Stxbp5l 
## Negative:  Nos1, Auts2, Etv1, Egfem1, Cadps2, Kcnq5, Gfra1, Schip1, Asic2, Dgkb 
##     Cmah, Dach1, Kcnt2, Epha5, Kcnab1, Rgs6, Stxbp6, Alcam, Creb5, Cntnap5a 
##     Il1rapl1, Lrrc4c, Hs6st3, Tmem108, Adgrl3, Cdh11, Gria3, Ebf1, Gramd1b, Adcy2 
## PC_ 3 
## Positive:  Cdh18, Klhl1, Pbx3, Kcnip4, Meis1, Sema5a, Kctd8, C79798, Gabrg3, Htr4 
##     Vwc2l, Gpc6, Cntn3, Dlc1, Serpini1, Zbbx, Csmd3, Zfhx3, Cdh9, March1 
##     Galnt18, Skap1, Plcl1, Stxbp5l, Pakap, Pde4d, Khdrbs2, Csmd1, Cadm2, Nrp2 
## Negative:  Adgrg6, Sgcd, Cysltr2, Slc4a4, Itgb6, Nfia, Nmu, Nos1, Gpr149, Grin3a 
##     Dapk2, Ano2, Rgs6, Cbln2, Lhfpl2, Zfp804a, Ccbe1, Dgkg, Zfp536, Kcnab1 
##     Otof, Efr3a, Cntnap5a, Cpne4, Zeb2, Smad6, Syt15, Gfra1, Scn11a, Zbtb7c 
## PC_ 4 
## Positive:  Dock10, Kcnip4, Lingo2, Csmd3, Sorcs1, Gda, Ndst4, Cntn5, Tac1, Nxph1 
##     Cntn3, Thsd4, Kctd8, Nyap2, Sgcz, Kcnq5, Lrp1b, Robo1, Plcb1, Pld5 
##     Tenm2, Lin7a, Lrrc7, Abca5, Lrrtm4, Dmd, Pgm5, Kctd16, Sorcs3, Prkg1 
## Negative:  Klhl1, Sema5a, Vwc2l, March1, Rasgef1b, Cdh9, Alk, Prune2, Galnt18, Il1rapl1 
##     Chsy3, Mgat4c, Kcnh7, Pbx3, Bcl2, Zfhx3, Dpp10, Zbbx, Galnt17, Grm5 
##     Pcdh7, C79798, Vcan, Kcnb2, Lncbate1, Auts2, Ghr, Dcc, Scgn, Olfr78 
## PC_ 5 
## Positive:  Prkg1, Kcnt2, Dgkb, Alcam, Car10, Rasgef1b, Galntl6, Epha5, Mgat4c, Vwc2l 
##     Vcan, Thsd7b, Cdh20, Dmd, Dpp10, Rgs6, Sema5a, Gucy1a2, Lingo2, Gabrg3 
##     Cdh9, Khdrbs2, Slc2a13, Kctd8, Lrrc4c, Htr4, Antxr2, Ptprz1, Adcy2, Man1a 
## Negative:  Ntng1, Ebf1, Trps1, Trhde, Cntn4, Csmd1, Gna14, Zmat4, Col18a1, Nkain3 
##     Tmtc2, Myo1e, Tenm4, Sctr, Nxph2, Cpa6, Nrg1, Kcnd2, Npas3, Neurod6 
##     Fstl4, Sez6l, Nav2, Ccdc85a, Myo16, Gal, Glp1r, Fbn2, Tfcp2l1, Nrp1

length(setdiff(VariableFeatures(object = GEX.seur),
                                                c(MT_gene,DIG, CC_gene) ))

## [1] 1067

head(setdiff(VariableFeatures(object = GEX.seur),
                                                c(MT_gene,DIG, CC_gene) ),300)

##   [1] "Mgat4c"     "Gal"        "Zfp804a"    "Grm5"       "Cntn5"     
##   [6] "Nmu"        "Klhl1"      "Adgrg6"     "Col24a1"    "Cpne4"     
##  [11] "Kcnb2"      "Robo2"      "Ntng1"      "Cntnap2"    "Tmeff2"    
##  [16] "Brinp3"     "Nrxn3"      "Lingo2"     "Ano2"       "Gpr149"    
##  [21] "Dgkg"       "Cdh8"       "Cdh18"      "Zfp804b"    "Ebf1"      
##  [26] "Cdh9"       "Ntrk3"      "Dgki"       "Sgcz"       "Pdzrn4"    
##  [31] "Cmah"       "Nxph2"      "Unc5d"      "Astn2"      "Kcnip4"    
##  [36] "Sema5a"     "Pcdh11x"    "Schip1"     "Pcdh9"      "Prkg1"     
##  [41] "Itgb6"      "Vwc2l"      "Cdh6"       "Nrg1"       "Postn"     
##  [46] "Dcc"        "Calcb"      "Car10"      "Piezo2"     "Grin3a"    
##  [51] "Myl1"       "March1"     "Clstn2"     "Vip"        "Rasgef1b"  
##  [56] "Bmpr1b"     "Efr3a"      "Ccbe1"      "Rab27b"     "Plxna4"    
##  [61] "Fut9"       "Dpp10"      "Gpc5"       "Pbx3"       "Asic2"     
##  [66] "Egfem1"     "Pcdh10"     "Abi3bp"     "Cysltr2"    "Serpini1"  
##  [71] "Col25a1"    "Skap1"      "Nkain3"     "Robo1"      "Fgf14"     
##  [76] "Trhde"      "Tafa2"      "Spock3"     "Chrna9"     "Zfhx3"     
##  [81] "Csmd3"      "Galnt13"    "Ghr"        "Tac1"       "Pde4d"     
##  [86] "Luzp2"      "Pde1a"      "Tbx20"      "Trps1"      "Kcnq5"     
##  [91] "Mid1"       "Cemip"      "Dgkb"       "Myo3b"      "Kcnh7"     
##  [96] "Acp7"       "Nos1"       "Flrt2"      "Plcl1"      "Iqgap2"    
## [101] "Cacna2d3"   "Csmd1"      "Tnr"        "Kctd16"     "Cadm2"     
## [106] "Parm1"      "Gna14"      "Cntn4"      "Chrm3"      "Nell1"     
## [111] "Abca9"      "Slc4a4"     "Il1rapl1"   "Nxph1"      "Kctd8"     
## [116] "Adamts6"    "Fbxl7"      "Myh3"       "Cartpt"     "Hs6st3"    
## [121] "Arhgap6"    "Kif26b"     "Penk"       "Aff2"       "Grm7"      
## [126] "Galnt18"    "Cacnb4"     "Synpr"      "Capn9"      "Sst"       
## [131] "Lama2"      "Ano5"       "Cdh10"      "Grik1"      "Opcml"     
## [136] "Kcnd2"      "Cadps2"     "Exoc3l2"    "Esr1"       "Galntl6"   
## [141] "L3mbtl4"    "Prune2"     "Gabrg3"     "Spag16"     "Epha7"     
## [146] "Cntn3"      "Prr16"      "Cux2"       "Calca"      "Chsy3"     
## [151] "Dach1"      "Antxr2"     "Morc1"      "Met"        "Thsd7b"    
## [156] "Kirrel3"    "Hs3st4"     "Zbbx"       "Npas3"      "Terb2"     
## [161] "Cbln2"      "Loxl1"      "Il18r1"     "Stab1"      "Galr1"     
## [166] "Runx1"      "Egfr"       "Ntrk2"      "Dock10"     "Olfr78"    
## [171] "Rasgrf1"    "Scnn1a"     "Scn7a"      "Cpa6"       "Col8a1"    
## [176] "Shisa9"     "D5Ertd615e" "C79798"     "Pde7b"      "Rerg"      
## [181] "Terf1"      "Ttc29"      "Pak7"       "Vcan"       "Etv1"      
## [186] "Grp"        "Olfr889"    "Adam2"      "Rarb"       "Tenm2"     
## [191] "Lhfpl2"     "Sorcs3"     "Gabrb1"     "Alcam"      "Bpifc"     
## [196] "Otof"       "Necab1"     "Gng2"       "Calcrl"     "Rxfp2"     
## [201] "Dapk2"      "Slco1a1"    "Ntm"        "Pdgfra"     "Plxdc2"    
## [206] "Lrrc4c"     "Fstl4"      "Bche"       "Mrc1"       "Nell1os"   
## [211] "Stxbp6"     "Hcn1"       "F13a1"      "Htr4"       "Nr4a3"     
## [216] "Elf5"       "Slc7a15"    "Myo1h"      "Fzd7"       "Thsd4"     
## [221] "Creb5"      "Sorcs1"     "Lockd"      "Tenm4"      "Tafa1"     
## [226] "Dysf"       "Rxfp1"      "Nwd2"       "Igfbp5"     "Rmst"      
## [231] "Gfra1"      "Dlc1"       "Ppp3ca"     "Plcb1"      "Rtl4"      
## [236] "Cfh"        "Cblc"       "Sctr"       "Arhgap15"   "Nfatc1"    
## [241] "Rgs6"       "Bmp5"       "Slc44a5"    "Ngfr"       "Ifi203"    
## [246] "Lamc3"      "Itgbl1"     "Nox4"       "Cntnap5b"   "Naaladl2"  
## [251] "Tex15"      "Htr2c"      "Adgrl2"     "Scube1"     "Rbpms"     
## [256] "Tex21"      "Casp8"      "Syt9"       "Etl4"       "Dab2"      
## [261] "Ptprt"      "Snx7"       "Tll1"       "Moxd1"      "Plac8"     
## [266] "Grid1"      "Agrn"       "Prkag2"     "Bcl2"       "Tenm1"     
## [271] "Prkca"      "C1qtnf7"    "Ror1"       "Dab1"       "Eya1"      
## [276] "Pde1c"      "Edil3"      "Hs3st2"     "Fam47e"     "Fbxw21"    
## [281] "Dgat2l6"    "Col18a1"    "Tmem64"     "Trpc3"      "Adam12"    
## [286] "Cfap299"    "Best2"      "Hsd17b14"   "Cpne8"      "Vit"       
## [291] "Adcy2"      "Myo1e"      "Cldn5"      "Rad51b"     "Brinp2"    
## [296] "Fras1"      "Fgf10"      "Zbtb7c"     "Grm8"       "Atpaf2"

DimPlot(GEX.seur, reduction = "pca",dims = 1:2, group.by = "orig.ident") +
  DimPlot(GEX.seur, reduction = "pca",dims = 3:4, group.by = "orig.ident")

DimPlot(GEX.seur, reduction = "pca",dims = 1:2, group.by = "FB.info", cols = color.FB) +
  DimPlot(GEX.seur, reduction = "pca",dims = 3:4, group.by = "FB.info", cols = color.FB)

DimHeatmap(GEX.seur, dims = 1:12, cells = 1500, balanced = TRUE,ncol = 4)

decide PCs to use

ElbowPlot(GEX.seur,ndims = 50)

PCs <- 1:18
GEX.seur <- FindNeighbors(GEX.seur, dims = PCs, k.param =20)

## Computing nearest neighbor graph

## Computing SNN

GEX.seur <- FindClusters(GEX.seur, method = 'igraph' ,resolution = 1.5)

## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 8342
## Number of edges: 333410
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8336
## Number of communities: 21
## Elapsed time: 0 seconds

Run UMAP/tSNE

GEX.seur <- RunTSNE(GEX.seur, dims=PCs, complexity=100)
GEX.seur <- RunUMAP(GEX.seur, dims=PCs, n.neighbors = 20, seed.use = 133)

## Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
## To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
## This message will be shown once per session

## 15:21:35 UMAP embedding parameters a = 0.9922 b = 1.112

## 15:21:35 Read 8342 rows and found 18 numeric columns

## 15:21:35 Using Annoy for neighbor search, n_neighbors = 20

## 15:21:35 Building Annoy index with metric = cosine, n_trees = 50

## 0%   10   20   30   40   50   60   70   80   90   100%

## [----|----|----|----|----|----|----|----|----|----|

## **************************************************|
## 15:21:36 Writing NN index file to temp file C:\Users\Shaorui\AppData\Local\Temp\Rtmp6hZbq2\file787045932716
## 15:21:36 Searching Annoy index using 1 thread, search_k = 2000
## 15:21:37 Annoy recall = 100%
## 15:21:38 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 20
## 15:21:39 Initializing from normalized Laplacian + noise (using irlba)
## 15:21:39 Commencing optimization for 500 epochs, with 239650 positive edges
## 15:22:00 Optimization finished

DimPlot(GEX.seur, reduction = "tsne", label = T) + DimPlot(GEX.seur, reduction = "umap", label = T)

DimPlot(GEX.seur, reduction = "umap", group.by = "preAnno1", label = T, label.size = 2.5) + 
  DimPlot(GEX.seur, reduction = "umap", group.by = "preAnno2", label = T, repel = F)

FeaturePlot(GEX.seur, reduction = "umap", features = c("nFeature_RNA","nCount_RNA","percent.mt","percent.rb"))

DimPlot(GEX.seur, label = F, pt.size = 0.05, repel = F, reduction = 'tsne', group.by = "FB.info", split.by = "FB.info", 
        ncol = 4, cols = color.FB)

DimPlot(GEX.seur, label = F, pt.size = 0.05, repel = F, reduction = 'umap', group.by = "FB.info", split.by = "FB.info", 
        ncol = 4, cols = color.FB)

refmapping

NatNeur2021

ref.seur <- readRDS("../analysis_ref/GSE149524.P21.integration_Anno.s.rds")
ref.seur

## An object of class Seurat 
## 37583 features across 4419 samples within 3 assays 
## Active assay: SCT (16365 features, 0 variable features)
##  2 other assays present: RNA, integrated
##  3 dimensional reductions calculated: pca, tsne, umap

ref.seur <- RunUMAP(ref.seur, assay = "integrated", umap.method = "uwot-learn", dims = 1:30)

## Warning in RunUMAP.default(object = data.use, reduction.model =
## reduction.model, : 'uwot-learn' is deprecated. Set umap.method = 'uwot' and
## return.model = TRUE

## UMAP will return its model

## 15:22:15 UMAP embedding parameters a = 0.9922 b = 1.112

## 15:22:15 Read 4419 rows and found 30 numeric columns

## 15:22:15 Using Annoy for neighbor search, n_neighbors = 30

## 15:22:15 Building Annoy index with metric = cosine, n_trees = 50

## 0%   10   20   30   40   50   60   70   80   90   100%

## [----|----|----|----|----|----|----|----|----|----|

## **************************************************|
## 15:22:15 Writing NN index file to temp file C:\Users\Shaorui\AppData\Local\Temp\Rtmp6hZbq2\file7870317045c9
## 15:22:15 Searching Annoy index using 1 thread, search_k = 3000
## 15:22:16 Annoy recall = 100%
## 15:22:17 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
## 15:22:17 Initializing from normalized Laplacian + noise (using irlba)
## 15:22:17 Commencing optimization for 500 epochs, with 181154 positive edges
## 15:22:30 Optimization finished

mapping

anchors

anchors.ref <- FindTransferAnchors(
  reference = ref.seur,
  reference.assay = "integrated",
  query = GEX.seur,
  #query.assay = "RNA",
  query.assay = "RNA",
  reference.reduction = 'pca',
  normalization.method = "LogNormalize",
  #normalization.method = "SCT",
  reduction = "pcaproject",
  dims = 1:30,
  k.anchor = 80,
  k.filter = 120,
  k.score = 60,
  max.features = 1600,
  nn.method = "annoy",
  n.trees = 100,
  eps = 0
)

## Projecting cell embeddings

## Finding neighborhoods

## Finding anchors

##  Found 20720 anchors

## Filtering anchors

##  Retained 6498 anchors

prediction

local_neur.pred <- MapQuery(
  anchorset = anchors.ref,
  query = GEX.seur,
  reference = ref.seur,
  refdata = list(
    ref.Anno1 = 'Anno1',
    ref.Anno2 = 'Anno2'
  ),
  reference.reduction = 'pca',
  reduction.model = 'umap'
)

## Warning: `invoke()` is deprecated as of rlang 0.4.0.
## Please use `exec()` or `inject()` instead.
## This warning is displayed once per session.

## Finding integration vectors

## Finding integration vector weights

## Predicting cell labels

## Warning: Keys should be one or more alphanumeric characters followed
## by an underscore, setting key from predictionscoreref.Anno1_ to
## predictionscorerefAnno1_

## Predicting cell labels

## Warning: Keys should be one or more alphanumeric characters followed
## by an underscore, setting key from predictionscoreref.Anno2_ to
## predictionscorerefAnno2_

## 
## Integrating dataset 2 with reference dataset

## Finding integration vectors

## Integrating data

## Computing nearest neighbors

## Running UMAP projection

## 15:24:14 Read 8342 rows

## 15:24:14 Processing block 1 of 1

## 15:24:14 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
## 15:24:14 Initializing by weighted average of neighbor coordinates using 1 thread
## 15:24:14 Commencing optimization for 167 epochs, with 250260 positive edges
## 15:24:22 Finished

local_neur.pred

## An object of class Seurat 
## 23284 features across 8342 samples within 3 assays 
## Active assay: RNA (23258 features, 1500 variable features)
##  2 other assays present: prediction.score.ref.Anno1, prediction.score.ref.Anno2
##  5 dimensional reductions calculated: pca, tsne, umap, ref.pca, ref.umap

check umap

#
raw1_p1 <- DimPlot(local_neur.pred,
                  reduction = 'umap',
                  group.by = 'preAnno1',
                  label=TRUE,
                  label.size = 3,
                  repel = F) + NoLegend()+ labs(title="Baf53Nb preAnno1")

raw1_p2 <- DimPlot(local_neur.pred,
                  reduction = 'umap',
                  group.by = 'predicted.ref.Anno1',
                  label=TRUE,
                  label.size = 3,
                  repel = T) + NoLegend()
raw1_p3 <- DimPlot(local_neur.pred,
                  reduction = 'umap',
                  group.by = 'preAnno2',
                  label=TRUE,
                  label.size = 3,
                  repel = TRUE) + NoLegend()+ labs(title="Baf53Nb preAnno2") 
raw1_p4 <- DimPlot(local_neur.pred,
                  reduction = 'umap',
                  group.by = 'predicted.ref.Anno2',
                  label=TRUE,
                  label.size = 3,
                  repel = TRUE) + NoLegend()

cowplot::plot_grid(raw1_p1,raw1_p2,ncol = 2)

cowplot::plot_grid(raw1_p3,raw1_p4,ncol = 2)

cowplot::plot_grid(FeaturePlot(local_neur.pred, features = "predicted.ref.Anno1.score", reduction = "umap"),
  FeaturePlot(local_neur.pred, features = "predicted.ref.Anno2.score", reduction = "umap"),ncol = 2)

newAnno

clusters sort new orders for mapping comparison,

and do newAnno based on it through ref.mapping to NatNeur2021,

newAnno would be like

EMN1-5 C0/1/3-C9-C2-C14-C16

IMN1-3 C5/4/15-C6-C11

IN1-3 C12-C18-C20(Sst+)

IPAN1-4 C7/10-C8/17-C19-C13

GEX.seur$sort_new <- factor(as.character(GEX.seur$seurat_clusters),
                            levels = c(0,1,3,
                                       9,
                                       2,
                                       14,
                                       16,
                                       5,4,15,
                                       6,
                                       11,
                                       12,
                                       18,
                                       20,
                                       7,10,
                                       8,17,
                                       19,
                                       13))

clusters to refAnno

#
local_neur.pred$sort_new <- GEX.seur$sort_new

pred.ref.Anno1 <- table(prediction=local_neur.pred$predicted.ref.Anno1,
                        clusters=local_neur.pred$sort_new)[paste0("ENC",c(1,#2,
                                                        3,4,
                                                        8,9,
                                                        10,#11,5,
                                                        6,7,12)),]
pred.ref.Anno1

##           clusters
## prediction   0   1   3   9   2  14  16   5   4  15   6  11  12  18  20   7  10
##      ENC1  943 884 683 288 125  12   4  63 120   0 185 145   0   1  24   1   0
##      ENC3    0   0   0 106 586 151   0   0   0   0   0   0   0   0   0   0   0
##      ENC4    0   0   0   0   0  55 200   0   0   0   0   0   0   0   0   0   0
##      ENC8    0   0   0   0   0   0   0 470 449 213 227   0   0   1   0   0   0
##      ENC9    0   0   0   0   0   0   0   0   1   0  63 123   0   0   0   0   0
##      ENC10   0   1   0   0  28  13   2  15  90   0  10  20 257 134   0   0   0
##      ENC6    0   0   0   0   0   0   0   0   0   0   0   0   0   0   1 423 316
##      ENC7    0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   5  11
##      ENC12   0   0   0   0   0   0   0   0   0   0   0   4   0   0   0   0   0
##           clusters
## prediction   8  17  19  13
##      ENC1    0   1   9   2
##      ENC3    0   0   0   0
##      ENC4    0   0   0   0
##      ENC8    0   0   1   0
##      ENC9    0   0   0   0
##      ENC10   0   0   7   3
##      ENC6    0   0   0   0
##      ENC7  405 178   1   0
##      ENC12   0   0  45 230

pred.ref.Anno2 <- table(prediction=local_neur.pred$predicted.ref.Anno2,
                        clusters=local_neur.pred$sort_new)
pred.ref.Anno2

##           clusters
## prediction   0   1   3   9   2  14  16   5   4  15   6  11  12  18  20   7  10
##      EMN1  943 884 683 290 197  25   4  63 120   0 185 145   0   1  24   1   0
##      EMN3    0   0   0 104 287  10   0   0   0   0   0   0   0   0   0   0   0
##      EMN4    0   0   0   0 222 111   0   0   0   0   0   0   0   0   0   0   0
##      EMN5    0   0   0   0   0  72 200   0   0   0   0   0   0   0   0   0   0
##      IMN1    0   0   0   0   0   0   0 470 449 213 227   0   0   1   0   0   0
##      IMN2    0   0   0   0   0   0   0   0   1   0  63 123   0   0   0   0   0
##      IN1     0   1   0   0  33  13   2  15  90   0  10  20 257 134   0   0   0
##      IN2     0   0   0   0   0   0   0   0   0   0   0   0   0   7   0   0   0
##      IPAN1   0   0   0   0   0   0   0   0   0   0   0   0   0   0   1 423 316
##      IPAN2   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   5  11
##      IPAN4   0   0   0   0   0   0   0   0   0   0   0   4   0   0   0   0   0
##           clusters
## prediction   8  17  19  13
##      EMN1    0   1   9   2
##      EMN3    0   0   0   0
##      EMN4    0   0   0   0
##      EMN5    0   0   0   0
##      IMN1    0   0   1   0
##      IMN2    0   0   0   0
##      IN1     0   0   7   3
##      IN2     0   0   0   0
##      IPAN1   0   0   0   0
##      IPAN2 405 178   1   0
##      IPAN4   0   0  45 230

cowplot::plot_grid(
pheatmap::pheatmap(pred.ref.Anno1,
                   main = "Cell Count (Baf53Nb mapping to NatNeur2021 P21, Anno1)",
      gaps_row = c(3,5,6),
      gaps_col = c(7,12,15),
                   cluster_rows = F, cluster_cols = F, display_numbers = T, number_format = "%.0f", legend = F, silent = T)$gtable,

pheatmap::pheatmap(t(100*rowRatio(t(pred.ref.Anno1))),
                   main = "Cell Ratio (Baf53Nb mapping to NatNeur2021 P21, Anno1)",
      gaps_row = c(3,5,6),
      gaps_col = c(7,12,15),
                   cluster_rows = F, cluster_cols = F, display_numbers = T, number_format = "%.1f", legend = F, silent = T)$gtable,
ncol = 1)

cowplot::plot_grid(
pheatmap::pheatmap(pred.ref.Anno2,
                   main = "Cell Count (Baf53Nb mapping to NatNeur2021 P21, Anno2)",
      gaps_row = c(4,6,8),
      gaps_col = c(7,12,15),
                   cluster_rows = F, cluster_cols = F, display_numbers = T, number_format = "%.0f", legend = F, silent = T)$gtable,

pheatmap::pheatmap(t(100*rowRatio(t(pred.ref.Anno2))),
                   main = "Cell Ratio (Baf53Nb mapping to NatNeur2021 P21, Anno2)",
      gaps_row = c(4,6,8),
      gaps_col = c(7,12,15),
                   cluster_rows = F, cluster_cols = F, display_numbers = T, number_format = "%.1f", legend = F, silent = T)$gtable,
ncol = 1)

newAnno

GEX.seur$newAnno <- as.character(GEX.seur$seurat_clusters)

GEX.seur$newAnno[GEX.seur$newAnno %in% c(0,1,3)] <- "EMN1"
GEX.seur$newAnno[GEX.seur$newAnno %in% c(9)] <- "EMN2"
GEX.seur$newAnno[GEX.seur$newAnno %in% c(2)] <- "EMN3"
GEX.seur$newAnno[GEX.seur$newAnno %in% c(14)] <- "EMN4"
GEX.seur$newAnno[GEX.seur$newAnno %in% c(16)] <- "EMN5"

GEX.seur$newAnno[GEX.seur$newAnno %in% c(5,4,15)] <- "IMN1"
GEX.seur$newAnno[GEX.seur$newAnno %in% c(6)] <- "IMN2"
GEX.seur$newAnno[GEX.seur$newAnno %in% c(11)] <- "IMN3"

GEX.seur$newAnno[GEX.seur$newAnno %in% c(12)] <- "IN1"
GEX.seur$newAnno[GEX.seur$newAnno %in% c(18)] <- "IN2"
GEX.seur$newAnno[GEX.seur$newAnno %in% c(20)] <- "IN3"

GEX.seur$newAnno[GEX.seur$newAnno %in% c(7,10)] <- "IPAN1"
GEX.seur$newAnno[GEX.seur$newAnno %in% c(8,17)] <- "IPAN2"
GEX.seur$newAnno[GEX.seur$newAnno %in% c(19)] <- "IPAN3"
GEX.seur$newAnno[GEX.seur$newAnno %in% c(13)] <- "IPAN4"

GEX.seur$newAnno <- factor(GEX.seur$newAnno,
                           levels = c(paste0("EMN",1:5),
                                      paste0("IMN",1:3),
                                      paste0("IN",1:3),
                                      paste0("IPAN",1:4)))

# define newAnno colors
color.newA <- c("#8AB6CE","#678BB1","#3975C1","#669FDF","#4CC1BD",
                "#BF7E6B","#D46B35","#FF8080",
                "#BDAE8D","#BD66C4","#C03778",
                "#97BA59","#DFAB16","#2BA956","#9FE727")

cowplot::plot_grid(
DimPlot(GEX.seur, reduction = "umap", group.by = "seurat_clusters", label = T, label.size = 3.25,repel = F, pt.size = 0.45),
  DimPlot(GEX.seur, label = F, pt.size = 0.45, repel = F, reduction = 'umap', group.by = "cnt",
        ncol = 1, cols = color.FB[c(5,1)]) ,
rel_widths = c(4.8,5),
ncol = 2)

cowplot::plot_grid(
DimPlot(GEX.seur, reduction = "umap", group.by = "newAnno", label = T, label.size = 3.25,repel = F, pt.size = 0.45),
  DimPlot(GEX.seur, label = F, pt.size = 0.45, repel = F, reduction = 'umap', group.by = "cnt",
        ncol = 1, cols = color.FB[c(5,1)]) ,
rel_widths = c(4.75,5),
ncol = 2)

cowplot::plot_grid(
DimPlot(GEX.seur, reduction = "umap", group.by = "newAnno", label = T, label.size = 3.25,repel = F, pt.size = 0.45,
        cols = color.newA),
  DimPlot(GEX.seur, label = F, pt.size = 0.45, repel = F, reduction = 'umap', group.by = "cnt",
        ncol = 1, cols = color.FB[c(5,1)]) ,
rel_widths = c(4.75,5),
ncol = 2)

Cell Composition

stat

GEX.seur@meta.data$cnt <- factor(gsub("1$|2$|3$|4$","",as.character(GEX.seur@meta.data$FB.info)))
GEX.seur$rep <- paste0("rep",
                        gsub("Baf53Nb.IF|Baf53Nb.CL","",as.character(GEX.seur$FB.info)))

stat_newAnno <- GEX.seur@meta.data[,c("cnt","rep","newAnno")]

stat_newAnno.s <- stat_newAnno %>%
  group_by(cnt,rep,newAnno) %>%
  summarise(count=n()) %>%
  mutate(prop= count/sum(count)) %>%
  ungroup

## `summarise()` has grouped output by 'cnt', 'rep'. You can override using the
## `.groups` argument.

pcom.newAnno <- stat_newAnno.s %>%
  ggplot(aes(x = newAnno, y = 100*prop, color=cnt)) +
  geom_bar(stat="summary", fun="mean", position = position_dodge(0.75), width = 0.58, fill="white") +
  theme_classic(base_size = 15) +
  scale_color_manual(values = c("skyblue","pink"), name="") +
  labs(title="Cell Composition of Baf53Nb data, newAnno", y = "Proportion") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6)) +
  
  geom_point(aes(x=newAnno, y = 100*prop, color= cnt),
             position = position_dodge(0.75),
             shape=16,alpha=0.75,size=2.15,
             stroke=0.15, show.legend=F)

pcom.newAnno

sig.test

glm.nb - anova.LRT

Sort.N <- c("Baf53Nb.CL","Baf53Nb.IF")

glm.nb_anova.newAnno <- list()

for(x1 in 1:2){
  for(x2 in 1:2){
    if(x2>x1){
      
      stat_newAnno.s_N <- stat_newAnno.s %>% filter(cnt %in% c(Sort.N[x1],Sort.N[x2]))
      
      total.N <- stat_newAnno.s_N %>%
        group_by(cnt, rep) %>%
        summarise(total=sum(count)) %>% ungroup() %>% as.data.frame()
      
      rownames(total.N) <- paste0(as.character(total.N$cnt),
                                  "_",
                                  as.character(total.N$rep))
      
      stat_newAnno.s_N$total <- total.N[paste0(stat_newAnno.s_N$cnt,
                                            "_",
                                            stat_newAnno.s_N$rep),"total"]
      
      glm.nb_anova.newAnno[[paste0(Sort.N[x1],"vs",Sort.N[x2])]] <- list()
      
      for(AA in levels(stat_newAnno.s_N$newAnno)){
        glm.nb_anova.newAnno[[paste0(Sort.N[x1],"vs",Sort.N[x2])]][[AA]] <- tryCatch({
          aaa <- stat_newAnno.s_N %>% filter(newAnno==AA)
          aaa.glmnb <- MASS::glm.nb(count ~ cnt + offset(log(total)), data = aaa, maxit=1000)
          aaa.anova <- anova(aaa.glmnb, test = "LRT")
          aaa.anova$'Pr(>Chi)'[2]
        }, error=function(e){
        
          tryCatch({
            aaa <- stat_newAnno.s_N %>% filter(newAnno==AA)
            aaa.glmnb <- MASS::glm.nb(count ~ cnt + offset(log(total)), data = aaa, maxit=100)
            aaa.anova <- anova(aaa.glmnb, test = "LRT")
            aaa.anova$'Pr(>Chi)'[2]
          }, error=function(e){
            
            tryCatch({
              aaa <- stat_newAnno.s_N %>% filter(newAnno==AA)
              aaa.glmnb <- MASS::glm.nb(count ~ cnt + offset(log(total)), data = aaa, maxit=10)
              aaa.anova <- anova(aaa.glmnb, test = "LRT")
              aaa.anova$'Pr(>Chi)'[2]
            }, error = function(e){
              NA
            })     
          })
        })
      }
      glm.nb_anova.newAnno[[paste0(Sort.N[x1],"vs",Sort.N[x2])]] <- unlist(glm.nb_anova.newAnno[[paste0(Sort.N[x1],"vs",Sort.N[x2])]]) 
    }
  }
}

glm.nb_anova.newAnno_df <- t(data.frame(Reduce(function(x,y){rbind(x,y)},glm.nb_anova.newAnno)))
rownames(glm.nb_anova.newAnno_df) <- names(glm.nb_anova.newAnno)
colnames(glm.nb_anova.newAnno_df) <- gsub("X","C",colnames(glm.nb_anova.newAnno_df))
glm.nb_anova.newAnno_df

##                               EMN1      EMN2     EMN3      EMN4      EMN5
## Baf53Nb.CLvsBaf53Nb.IF 0.008660417 0.3309766 0.891037 0.8971586 0.6479739
##                              IMN1      IMN2       IMN3       IN1       IN2
## Baf53Nb.CLvsBaf53Nb.IF 0.07575576 0.1147984 0.09201819 0.8945842 0.5844769
##                              IN3      IPAN1     IPAN2     IPAN3     IPAN4
## Baf53Nb.CLvsBaf53Nb.IF 0.8401056 0.09909844 0.4667117 0.9061248 0.7642453

markers

long

markers.new <- list(EMN=c("Chat","Gfra2","Casz1","Xylt1",
                       "Ptprt","Bnc2","Tox","Oprk1",
                       "Kalrn","Pi15","Drd2","Adamtsl1", 
                       "Fbxw15","Fbxw24","Cdc14a","Chrna7",
                       "Caln1","Tmem132c","Satb1","Cntnap5b",
                       "Gabrb1","Nxph1","Lama2","Lrrc7",
                       "Ryr3","Eda","Chgb","Pgm5",
                       "Shc4","Vgll3","Ptn","Tac1",
                       "Kctd8","Ntrk2","Penk","Pde7b",
                       "Fut9","Nfatc1","Egfr","Mgll",
                       "Ntn1","Prlr","Chrm3"
                       ),
                 IMN=c("Nos1","Dach1","Kcnab1",
                       "Gfra1","Etv1",
                       "Man1a","Airn","Enpp1","Unc13c",
                       "Plpp3","Fat1","Adcy2","S100a4",
                       "Npy","Cmah","Creb5","Vip","Pde1a",
                       "Ebf1","Gpc5","Col25a1","Cartpt",
                       "Lgr5","Gabrb2"
                       ),
                 IN=c("Npas3","Synpr","St18","Gal",
                      "Kcnk13",
                      "Neurod6","Moxd1","Sctr","Fndc1",
                      "Piezo1","Sst","Adamts9","Kcnn2",
                      "Pantr1","Vwc2","Vipr2","Tacr1",
                      "Calb2"),
                 IPAN=c("Calcb","Nmu","Kcnj12","Nog",
                        "Bmp4","Adgrg6","Cysltr2","Pcdh10",
                        "Ngfr","Galr1","Met",
                        "Htr3a","Il7","Aff2","Gpr149",
                        "Efr3a","Cdh6","Cdh8","Pdzrn4",
                        "Clstn2","Cachd1","Ano2","Ntrk3",
                        "Cpne4","Vwc2l","Cdh9","Scgn",
                        "Vcan","Cck","Kcnh7","Piezo2",
                        "Abca6","Fam107b","Npy1r","Abca9",
                        "Abca8b","Rerg","Bmpr1b","Skap1",
                        "L3mbtl4","Tafa2","Nxph2","Gm32647",
                        "Gm29521","Ntng1"))


pm.CL_new <- DotPlot(subset(GEX.seur,subset=cnt %in% c("Baf53Nb.CL")), features = as.vector(unlist(markers.new)), group.by = "newAnno",
        cols = c("midnightblue","darkorange1")) +
  theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1,size = 9.15))+ scale_y_discrete(limits=rev) + labs(title="CL only")
pm.CL_new

pm.All_new <- DotPlot(GEX.seur, features = as.vector(unlist(markers.new)), group.by = "newAnno",
        cols = c("midnightblue","darkorange1")) +
  theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1,size = 9.15))+ scale_y_discrete(limits=rev) + labs(title="All")
pm.All_new

short

markers.new.s <- list(EMN=c("Chat","Bnc2","Tox","Ptprt",
                       "Gfra2","Oprk1","Adamtsl1", 
                       "Fbxw15","Fbxw24","Chrna7",
                       "Satb1","Cntnap5b",
                       "Gabrb1","Nxph1","Lama2","Lrrc7",
                       "Ryr3","Eda","Tac1",
                       "Kctd8","Ntrk2","Penk",
                       "Fut9","Nfatc1","Egfr","Mgll",
                       "Chrm3"
                       ),
                 IMN=c("Nos1","Kcnab1",
                       "Gfra1","Etv1",
                       "Man1a","Airn",
                       "Adcy2","Cmah","Creb5","Vip","Pde1a",
                       "Ebf1","Gpc5"
                       ),
                 IN=c("Npas3","Synpr","St18","Gal",
                      "Kcnk13",
                      "Neurod6","Moxd1","Sctr",
                      "Piezo1","Sst","Adamts9","Kcnn2",
                      "Calb2"),
                 IPAN=c("Calcb","Nmu","Adgrg6","Pcdh10",
                        "Ngfr","Galr1","Il7","Aff2",
                        "Gpr149","Cdh6","Cdh8",
                        "Clstn2","Ano2","Ntrk3",
                        "Cpne4","Vwc2l","Cdh9","Scgn",
                        "Vcan","Cck","Piezo2","Kcnh7",
                        "Rerg","Bmpr1b","Skap1","Ntng1",
                        "Tafa2","Nxph2"))


pm.CL_new.s <- DotPlot(subset(GEX.seur,subset=cnt %in% c("Baf53Nb.CL")), features = as.vector(unlist(markers.new.s)), group.by = "newAnno",
        cols = c("midnightblue","darkorange1")) +
  theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1,size = 9.15))+ scale_y_discrete(limits=rev) + labs(title="CL only")
pm.CL_new.s

pm.All_new.s <- DotPlot(GEX.seur, features = as.vector(unlist(markers.new.s)), group.by = "newAnno",
        cols = c("midnightblue","darkorange1")) +
  theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1,size = 9.15))+ scale_y_discrete(limits=rev) + labs(title="All")
pm.All_new.s

#saveRDS(GEX.seur,"./Baf53Nb.Ileum_IFd5.pure_newAnno20230728.rds")